arXiv:1502.01260v2 [stat.ME] 20 Oct 2015 


1 


Hyperspectral Unmixing with Spectral Variability 
Using a Perturbed Linear Mixing Model 

Pierre-Antoine Thouvenin, Student Member, IEEE, Nicolas Dobigeon, Senior Member, IEEE and 

Jean-Yves Tourneret, Senior Member, IEEE 


Abstract —Given a mixed hyperspectral data set, linear un¬ 
mixing aims at estimating the reference spectral signatures 
composing the data - referred to as endmembers - their abun¬ 
dance fractions and their number. In practice, the identified 
endmembers can vary spectrally within a given image and can 
thus be construed as variable instances of reference endmem¬ 
bers. Ignoring this variability induces estimation errors that 
are propagated into the unmixing procedure. To address this 
issue, endmember variability estimation consists of estimating 
the reference spectral signatures from which the estimated 
endmembers have been derived as well as their variability with 
respect to these references. This paper introduces a new linear 
mixing model that explicitly accounts for spatial and spectral 
endmember variabilities. The parameters of this model can 
be estimated using an optimization algorithm based on the 
alternating direction method of multipliers. The performance of 
the proposed unmixing method is evaluated on synthetic and real 
data. A comparison with state-of-the-art algorithms designed to 
model and estimate endmember variability allows the interest of 
the proposed unmixing solution to be appreciated. 

Index Terms —Hyperspectral imagery, linear unmixing, end- 
member spatial and spectral variability. Alternating Direction 
Method of Multipliers (ADMM). 

1. Introduction 

VER the past decades, hyperspectral imagery has been 
receiving an increasing interest. Whereas traditional 
red / green / blue or multispectral images are composed of 
a limited number of spectral channels (from three to tens), 
hyperspectral images are acquired in hundreds of contiguous 
spectral bands facilitating the analysis of the elements in the 
scene, e.g., determining their nature and relative proportions. 
However, the high spectral resolution of these images is 
mitigated by their lower spatial resolution, hence the need to 
unmix the data. Spectral unmixing is aimed at estimating the 
reference spectral signatures - referred to as endmembers - 
their abundance fractions and their number from which the L- 
multi-band observations are derived according to a predefined 
mixing model. Assuming the absence of any microscopic 
interaction between the materials of the imaged scene and a 
negligible declivity, a linear mixing model (LMM) is clas¬ 
sically used to describe the structure of the collected data 
[1]. However, the spectral signatures contained in a reference 
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image can vary spectrally, spatially or temporally from an 
image to another due to varying acquisition conditions. This 
can result in significant estimation errors being propagated 
throughout the unmixing process. Various models either de¬ 
rived from a statistical or a deterministic point of view have 
been designed to address this issue [2]. More precisely, the 
first class of methods assumes that the endmember spectra 
can be considered as realizations of multivariate distributions. 
The most popular models are the normal composition model 
[3] and the beta compositional model [4]. The second class 
of methods considers the endmember signatures as members 
of spectral libraries associated with each material (bundles). 
Two methods using spectral libraries have been especially 
considered in the literature: the automated endmember bundles 
(AEB) [5] and the Eisher discriminant null space (EDNS) 
[6]. Whereas AEB enables the extraction of an endmember 
library to account for spectral variabilities, the aim of EDNS 
is to estimate a transformation projection matrix to project 
the hyperspectral data into a space minimizing the variability 
impact. 

Since the identified endmembers can be considered as 
variable instances of reference endmembers, we introduce an 
extended version of the classical LMM to explicitly model 
the spectral variability. In [7], the variability is assumed 
to only result from scaling factors. Conversely, in this pa¬ 
per, inspired by a model designed in [8], each endmember 
is represented by a "pure" spectral signature corrupted by 
an additive perturbation accounting for its variability. The 
perturbation is allowed to vary from a pixel to another to 
represent spatial-spectral variabilities. As a result, the proposed 
perturbed LMM (PLMM) can capture endmember spatial and 
spectral variability within a given image. To the best of our 
knowledge, it is the first time endmember variability has been 
explicitly modeled as an additive perturbation. 

The promising results obtained with the alternating direction 
method of multipliers (ADMM) in hyperspectral imagery [9] 
and in image deblurring [10]-[14] serve as an incentive to 
apply a similar framework to conduct PLMM-based unmixing. 
A key property of the ADMM framework lies in the introduc¬ 
tion of appropriate splitting variables. Indeed, the specified 
constraints can be handled independently from the rest of the 
problem and often lead to analytical solutions when solving the 
resulting optimization problem. Using this fruitful principle, 
an ADMM-based algorithm for linear unmixing using a group 
lasso ^ 2 , 1 -norBi regularization was recently developed in [15], 
[16]. Inspired by these examples, this paper proposes to 
exploit the advantages of an ADMM-based resolution of the 
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TABLE I 

Table of notations. 


N 

number of pixels 

L 

number of spectral bands 

K 

number of endmembers 

Y e 

lexicographically ordered pixels 

M e 

endmember matrix 

dMn e 

nth variability matrix 

A gRifxAT 

abundance matrix 

V G 

projector on the space spanned by the 
K — 1 principal components of Y 

U G 

inverse projector 

TeR(^f-l)xif 

projection of M on the PC A space 


column k, row j of the matrix Mn 

column k, row j of the matrix AB 


term-wise inequality 

C+ 

^n,p 

= {X G R^-Xpjx h 0„,p} 

In 

= [1,1,...,!]^ GK'* 

_ f 0 ifx G A 

Xa{x) 

~ \ +00 otherwise. 


linear unmixing problem to account for spatial and spectral 
endmember variabilities. 

Throughout the article, the number of endmembers will be 
assumed to be a priori known or estimated by any state-of-the- 
art method (e.g., [17]) since estimating the required number 
of endmembers to appropriately describe the data as well as 
endmember variability is a challenging task. Indeed, the choice 
of K drastically alters the representation of the imaged scene, 
and is thus a crucial step to the endmember identification and 
the subsequent abundance estimation [1], [17], [18]. 

The paper is organized as follows. The PLMM accounting 
for spectral and spatial variabilities is introduced in Section II. 
Section III describes an ADMM-based algorithm to solve the 
resulting optimization problem. Experimental results obtained 
on synthetic and real data are reported in Section IV and V 
respectively. The results obtained with the proposed algorithm 
are systematically compared to those of the vertex component 
analysis / fully constrained least squares (VCA/FCLS), the 
simplex identification via split augmented Lagrangian (SISAL) 
[9] coupled with FCLS, AEB and FDNS. Section VI finally 
concludes this work. 

II. Problem statement 
A. Perturbed linear mixing model (PLMM) 

In the absence of any specific prior knowledge on the 
variability nature (i.e., errors affecting the endmembers), we 
have chosen to explicitly represent the variability by a spa¬ 
tially varying additive endmember perturbation. This choice, 
inspired by a model designed in [8], appears to be simple 
and fiexible enough to account for the observed variability. 
Assuming that the number of endmembers K is known, the 
proposed PLMM differs from the classical EMM insofar as 
each pixel is represented by a combination of the K 
endmembers - denoted as m/c - affected by a perturbation 
vector dnin k accounting for endmember variability. The 
resulting PLMM can be written 

K 

y„ = ^ akn (mfe + + b„ for n = 1,..., (1) 

fe=i 


where denotes the nth image pixel, mk is the /cth end- 
member, Okn is the proportion of the k\h endmember in the 
nth pixel, and dm^denotes the perturbation of the kih 
endmember in the nth pixel. Finally, models a zero-mean 
white Gaussian noise resulting from the data acquisition as 
well as modeling errors. We can note that the proposed PLMM 
reduces to the classical LMM in absence of variability. In 
matrix form, the PLMM (1) can be written as follows 


Y = MA- 


dMiai 




+B (2) 


where Y = [yi,...,yA] is an L x N matrix containing 
the image pixels, M is an L x AT matrix containing the 
endmembers, A is a A x Y matrix composed of the abundance 
vectors a^, dM^ is an L x A matrix whose columns are the 
perturbation vectors associated with the nth pixel, and B is an 
Lx N matrix accounting for the noise. The non-negativity and 
sum-to-one constraints usually considered to refiect physical 
considerations are 

A ^ Ok,n, = 1a 

M ^ 0l,a 5 M + dM^ ^ 0 l,A 5 Vn = 1,..., Y. 


When compared to the underlying models proposed in the 
literature to mitigate variability [2], model (1) presents the 
advantage to explicitly address the variability phenomenon in 
terms of an additive perturbation affecting each endmember. 
This perturbation accounts for any deviation from the linear 
mixing model (as will be illustrated in our experiments). The 
main contribution of this paper is to propose an unsupervised 
algorithm for estimating the endmembers contained in the 
image and the abundances and endmember variability for each 
pixel of this image. 


B. Problem formulation 

As mentioned in Section I, the PLMM (1) and constraints 
(3) can be combined to formulate a constrained optimization 
problem. An appropriate cost function is required to estimate 
the parameters M, A, dM. Assuming the signal is corrupted 
by a zero-mean white Gaussian noise, we define the data fitting 
term as the Frobenius norm of the difference between the 
acquisitions Y and the reconstructed data MA + A. Since the 
problem is ill-posed, additional penalization terms are needed. 
In this paper, we propose to define penalization functions 
and T to refiect the available a priori knowledge on M, A 
and dM respectively. As a result, the optimization problem is 
expressed as 

(M*,dM*, A*) G arg min |j(M,dM, A) s.t. (3)| (4) 

M,dM,A ^ ^ 

with 

J(M, dM, A) ||Y - MA - A||2 + a^A)+ 

where the penalization parameters control the trade¬ 

off between the data fitting term ^ ||Y — MA — A||p and the 
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penalties ^(A), ^(M) and T(dM). In addition, we assume 
that the penalization functions are separable, leading to 

N 

i>(A) = ^ <^(a„) (6) 

n=l 

L 

= (7) 

i=l 

N 

T(dM) = ^^;(dM„) (8) 

n=l 

where denotes the ^th row of M and 0, t/; and v are 
non-negative differentiable convex functions. This assumption 
is used to decompose (4) into a collection of simpler sub¬ 
problems described in Section III. All these penalizations are 
described in the next paragraphs. 

1) Abundance penalization: The abundance penalization ^ 
has been chosen to promote spatially smooth abundances as 
in [19]. More precisely, the abundance spatial smoothness 
penalization is expressed in matrix form as 

$(A) = \ ||AH||^ (9) 

where H is a matrix computing the differences between 
the abundances of a given pixel and those of its 4 nearest 
neighbors [19]. The resulting expression of (j) is detailed in 
Appendix A. 

2) Endmember penalization: As for classical penaliza¬ 
tions found in the literature consist of constraining the size 
of the simplex whose vertices are the endmember signatures. 
The volume criterion used in [20], [21] enables the volume 
exactly occupied by the {K — 1)-simplex formed by the 
endmembers to be penalized. The mutual distance between 
the endmembers introduced in [22], [23] (which approximates 
the volume) has a similar purpose. Finally, if the endmembers 
are a priori close from available reference spectral signatures, 
a penalization on the distance between the endmembers and 
these signatures can be implemented. The expression of the 
distance between the endmembers and some reference spectral 
signatures, the mutual distance between the endmembers and 
the volume penalization are recalled in the following lines. For 
each penalization type, the corresponding expression of is 
given in Appendix A. 

• The distance between the endmembers and some refer¬ 
ence spectral signatures Mq is given by 

^'(M) = i||M-Mo||^ (10) 

• The mutual distance between the endmembers is ex¬ 
pressed in matrix form as 

i = l \j = l 

• Under the pure pixel and linear mixture assumptions, the 
data points are enclosed in a (A — 1)-simplex whose 
vertices are the endmembers [21]. Let T be the projection 
of M on the space spanned by the A — 1 principal 



Algorithm 1: PLMM-unmixing: global algorithm. 


Data: Y, A^o),M^o),dM(°) 

Result: A,M,dM 
begin 
k 1\ 

while stopping criterion not satisfied do 

AW ^ arg min a) ; 


MW ^ arg min j(^M,A^) ; 


^ arg min j(M('=),dM,AW) ; 


A: fc + 1; 

A^ AW; 

M ^ mW; 
dM ^ dM^*^ 


components of Y. The expression of the volume of this 
subspace is 


V(T) = 


1 

( 7 ^- 1 )! 



To ensure the differentiability of the penalization with 
respect to T, we propose to consider the following 
penalty 

,2(rr\ ( 12 ) 


1 , 


4f(M) = -V%T). 


3) Variability penalization: The variability penalizing func¬ 
tion T has been designed to limit the norm of the spectral 
variability. Indeed, it is legitimate to penalize the energy of 
the perturbation matrices dM^ in order to obtain a reasonable 
endmember variability. In this paper, we propose to consider 
the following penalty (having the advantage to be differen¬ 
tiable with respect to dM^) 


1 1 

T(dM) = -||dM||^ = -^||dM„||^ (13) 

n=l 

To the best of our knowledge, no specific information re¬ 
garding the spatial distribution of the variability is available 
in the remote sensing literature so far. We have consequently 
preferred not to include any additional regularization on dM. 
However, any spatial penalization satisfying the assumptions 
given in Paragraph II-B can be added when necessary (e.g., 
a group-Lasso ^ 2,1 penalization to promote spatial sparsity of 
the variability term dM). 


III. An ADMM-based algorithm 

Since the problem (4) is not convex, a minimization strategy 
similar to [12] has been adopted. Precisely, the cost function J 
is successively minimized with respect to each variable A, M 
and dM until a stopping criterion is satisfied. The assumptions 
made on the penalization functions 4>, T in Section II allow 
the global optimization problem to be divided into a collection 
of strictly convex sub-problems. These sub-problems have 
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the nice property to involve differentiable functions simpli¬ 
fying their resolution. Having introduced appropriate splitting 
variables to account for the constraints, these sub-problems 
are finally solved using ADMM steps admitting closed-form 
expressions due to the separability assumption. The three 
minimization steps considered in this algorithm present a 
highly similar structure. The details are reported in Appendix 
B to facilitate the reading of this paper. 


A. ADMM: general principle 

The ADMM is a technique combining the benefits of aug¬ 
mented Lagrangian and dual decomposition methods to solve 
constrained optimization problems [24]. More precisely, the 
method consists of solving the original optimization problem 
by successively minimizing the cost function of interest with 
respect to each variable. The following elements (extracted 
from [24]) recall a general formulation of the problem. Given 
/ : ^ M+, G ^ M+, A G and B G 

consider the general optimization problem 


Algorithm 2: ADMM optimization w.r.t. A (step (a)). 

Data; Y, eauai, 

Result: A 

for n = 1 to A do 

k ^ 1; 

= 0 ; 

while stopping criterion not satisfied do 
arg min £^^A)(fe-i) 

arg min £ (ax^-d ; 

w(A) \ ) 


(A)(fc) 


^ using (20) ; 

A: <(- fc + 1 ; 


- s ; 


, (k) 

a„ <(- a)i 


min 

x,z 


/(x) + g{z) 


Ax + Bz = c 


(14) 


The scaled augmented Lagrangian associated with this prob¬ 
lem can be written 


Cp (x, z, u) 

= /(x) + g{z) + ^ ||Ax + Bz - c + u ||2 

where p > 0. Denote as and the primal 

variables and the dual variable at iteration /c + 1 of the 
algorithm 


e arg min £p (x,z('=),u('')) 

^(k+i) 

e arg min £p (x(''+i),z,u('')) 

u(^+l) 

= + Ax(''+^) + Bz('=+^) - c. 


The ADMM consists in successively minimizing Cp with 
respect to x, z and u. A classical stopping criterion involves 
the primal and dual residuals at iteration k 1 (see [24, p. 
19]): the procedure is iterated until 


r 


(k) 


< and 



< 


(15) 


2 2 

where the primal and dual residuals at iteration k 1 are 
respectively given by 


j.ik+ 1 ) ^ - c (16) 

s(*+i) = pA'^B - z^'')) (17) 


and 


£Pri = ^rel ^ ^ ||^,||2| 


gdual ^ ^^abs ^rel 


A^y(^) 


(18) 

(19) 


Finally, the parameter p can be adjusted using the rule de¬ 
scribed in [24, p. 20] 




= < 


^incyfe) if||rW|| >p||s 


p 


p(fe)/^decr if ||sW ||2 > p||r 
Otherwise. 


W| 

ik)\ 


( 20 ) 


Note that this parameter adjustment does not alter the ADMM 
convergence as long as it is performed finitely many times. 


B. Optimization with respect to A 


With the assumptions made in paragraph II-B, optimizing 
the cost function J with respect to A under the constraints 
(3) is equivalent to solving the following problems 


a* = arg min 



(M + dM^)an||2 + 
an ^ Ok, aJlK = 1 


( 21 ) 

After introducing the splitting variables G R^ for 

n = 1,..., A such that 



( 22 ) 


the resulting scaled augmented Lagrangian is expressed as 


C 


4A) i ||y^ 


l^ 


(A) 


2 

■ Q(0(an) 


Qa„ + — s + A^^^ 


'^s+ 


.(w-«) 


{M + dM„)a„||, 

2 
2 


(23) 


with /in >0. The resulting algorithm (step a of Algo. 1) is 
detailed in Algo. 2, and the solution to each sub-problem is 
given in Appendix B. 
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C. Optimization with respect to M 


E. Convergence and computational cost 


Similarly to Paragraph III-B, optimizing J with respect to 
M under the constraint (3) is equivalent to solving 

2 


= arg mm < 


yi - meA - d. 
s.t. for n = 1 


+ Pipime) 


,N 


^ 0^, + dvcvn,i h 0 


K 


(24) 

where denotes the ^th row of M. Introducing the splitting 
variables G for ^ = 1,..., L such that 


1 

Iat 


fi., - wi“’ = - 


0l 


. \dmAr,<'/ 


(25) 


the associated scaled augmented Lagrangian can be written 

2 


jC (m) 






erne — W 


(M) 


yi - m^A - Si 
2 
F 


■Ar 


(26) 




with p,f^^ > 0. The resulting algorithm (step b of Algo. 1) is 
similar to Algo. 2. The solution to the optimization problems 
depends on the selected endmember penalizing function 
chosen in paragraph II-B2 (see Appendix B for more details). 

D. Optimization with respect to dM 

Finally, optimizing J with respect to dM under the con¬ 
straint (3) is equivalent to solving the sub-problems 
f i||y„-(M + dM„)a„||2 

+7r;(dM„) ^ . (27) 

S.t. M + dM^ ^ 

Introducing the splitting variables = M + dM^ for 

n = 1,..., A, the resulting scaled augmented Lagrangian is 
given by 


dM* = arg min < 

dM^ 


+ 


,,(dM) 

pn 


dM^ + M - + aI‘1’^) 


■ 7u(dM„) 




(28) 

with > 0. The resulting algorithm (step c of Algo. 1) 

is similar to Algo. 2. The solution to these problems is given 
in Appendix B. 


The optimization procedures detailed above are performed 
sequentially until the stopping criterion is satisfied. The next 
sections evaluate the performance of the resulting unmixing 
strategy via several experiments conducted on synthetic and 
real data. 


The alternating scheme proposed in Alg. 1 is nothing but a 
block coordinate descent descent (BCD) which is guaranteed 
to converge to a stationary point of the objective function 
J as long as each sub-problem is exactly minimized [25, 
Proposition 2.7.1]. Besides, the sub-problems tackled in III-B, 
III-C and III-D are strongly convex, hence the convergence 
of the ADMM steps toward the unique minimum of each 
independent sub-problem when the augmented Lagrangian 
parameter has a constant value (see for instance [24]). The 
same convergence result applies to the ADMM with the 
parameter adjustment introduced in Paragraph III-A as long 
as the parameter is updated finitely many times [24]. We 
may however mention that the proximal alternating linearized 
minimization (PALM) [26] could also be directly applied to the 
considered problem with a rigorous convergence proof based 
on the Kurdyka-Lojasiewicz property. This alternative work 
has been presented in [27]. 

Considering the significant number of unknown parameters 
and the simple expression of the ADMM updates detailed 
in Appendix B, we can note that the computational cost is 
dominated by matrix products, yielding an overall 0{LK‘^N) 
computational cost. 

IV. Experiment with synthetic data 

This section considers four images of size 128 x 64 
acquired in 413 bands. Each image corresponds to a mixture 
of K endmembers with K G {3,6} in presence or absence 
of pure pixels (the absence of pure pixels is considered to 
evaluate the algorithm performance in a very challenging 
scenario). The synthetic linear mixtures have been corrupted 
by additive white Gaussian noise to ensure the signal-to-noise 
ratio is SNR = 30dB. Since no accepted variability model is 
available in the literature, we propose the following generation 
procedure to introduce controlled spectral variability. The cor¬ 
rupted endmembers involved in the mixture (see Fig. 2) have 
been generated using the product of reference endmembers 
with randomly drawn piece-wise affine functions, providing 
realistic perturbed endmembers as represented in Fig. 1. For 
a given variability coefficient Cyar > 0 fixed by the user, the 
parameters ^i,i G {1,2,3} and Lbreak ^ {1? • • • ^ ^} introduced 
in Fig. 1 have been generated as follows 

^ ^l-c,^/2,l + Cvar/2] : ^ ^ {1: 2, 3} (29) 

inbreak = L^/2 + [LU/3\ J , f/ - A(0, 1) (30) 

where [-J denotes the floor operator. The synthetic data used 
in the proposed experiments have been generated with a 
value of Cyar that is lower in the upper half of the image 
(Cyar = 0.1) than in the lower half (Cyar = 0.25). Some 
instances of the corresponding perturbed endmember spectra 
are depicted in Fig. 2. Note that different affine functions have 
been considered for different endmembers and different pixels. 

A. State-of-the-art methods 

The results of the proposed algorithm have been compared 
to those obtained with two classical linear unmixing methods 
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Fig. 1. Example of a randomly-generated affine function used to generate 
the synthetically perturbed endmembers. 


(VCA/FCLS, SISAL/FCLS) and two variability accounting 
algorithms (AEB, FDNS). These methods are recalled below 






Classical unmixing methods (no variability) 

0.8 


/ 

0.4 


• VCA/FCLS: the endmembers are first extracted 
using the vertex component analysis [28]. The abun¬ 
dances are then estimated for each pixel using the 
fully constrained least squares (FCLS) algorithm 

OJ 

8 0.6 
ro 

8 

*§0.4 

0.2 

aJ 

kJ I 

8 

|0.3 

<u 

°^0.2 

0.1 

f \ 

[29]); 

• SISAL/FCLS: the endmembers are first extracted 

500 

1000 1500 2000 

X (nm) 

500 1000 1500 2000 

X (nm) 


using the simplex identification via split augmented 
Lagrangian [9]. The tolerance for the stopping rule 
has been set to 10“^ and VGA has been used 
as an initialization step. The abundances are then 
estimated for each pixel using FCLS. 


Fig. 2. Reference endmembers (red lines) and 20 corresponding instances 
under spectral variability (blue dotted lines) involved in the synthetic data 
experiments. 

TABLE II 

ADMM PARAMETERS. 


2) Variability accounting unmixing methods 

• AEB [5], [30], [31]: the size of the bundles is equal 
to 25% of the total pixel number. The endmembers 
and abundance are estimated using VCA/FCLS; 

• FDNS [6]: the endmembers and abundances are 
estimated by VCA/FCLS ; 

• Proposed method (BCD/ADMM): endmembers and 
abundances have been initialized with VCA/FCLS 
estimates. Note that VCA/FCLS is a method assum¬ 
ing that there are pure pixels in the image, which 
can be problematic in case the imaged scene does 
not satisfy this assumption. The variability matrices 
have been initialized with all their entries equal to 
eps^. The algorithm is stopped when the relative 
difference between two successive values of the 
objective function is less than 10“^. 

Different penalization combinations have been compared 
for the proposed method. The abbreviations ss, mv and 
vca are used for spatially smooth, minimum volume and 
minimum distance to VCA in the following. The absence of 
any additional abbreviation means that the method does not 
include any abundance or endmember penalization term. 



Synthetic data 

Real data 

.^incr 

1.1 

1.1 

^decr 

1.1 

1.1 

n 

10 

10 

(A)(0) 

10-4 

10-"^ 

,/M)(0) 

Up 

tdM)(0) 

Un 

10-® 

10-4 

10“® 

10-4 

^abs 

10-1 

(M 

1 

O 

1—1 

^rel 

1 

o 

1—1 

1 

o 

1—1 


angle mapper (aSAM) 

if ^ ||mfe||2||mfc||2 

as well as in terms of abundance and perturbation estimation 
by global mean square errors (GMSEs) 

1 


GMSE(A) = 
GMSE(dM) = 


KN 

1 

~ni7k 


lA-AII 


N 


^||dM„-dM, 


nllF- 


n=l 


The performance of the algorithm has been assessed in 
terms of endmember estimation using the average spectral 

^Matlab constant eps = 2.22 x 10“^®. 


As a measure of fit, the following reconstruction error (RE) 
has been also considered 



Y-Y 


2 

F 
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(a) Th. Abundance 1 (b) Th. Abundance 2 (c) Th. Abundance 3 



(d) Abundance 1 (e) Abundance 2 (f) Abundance 3 



(g) Variability 1 (h) Variability 2 (i) Variability 3 


Fig. 3. True abundances (fig. 3a to 3c) and the ssmvBCD/ADMM-estimations 
(fig. 3d to 4c) - obtained with a synthetic dataset (no pure pixels, K = 3). 
The spatial distribution of the variability with respect to each endmember is 
presented in terms of energy (^^ ||dm^^/e ||2 for the kth endmember in the 
nth pixel) for visualization purpose in Figs. 3g to 3i. 


where Y is the matrix formed of the pixels reconstructed using 
the parameters estimated by the algorithm. 


B. Results 

The parameters used for the ADMM algorithms are detailed 
in Table II, and the values chosen by cross-validation for a, 
P and 7 are reported in Table III and IV. The performance 
measures returned by the unmixing methods are provided in 
Table III for the datasets containing pure pixels, and in Table 
IV for images without pure pixels, leading to the following 
conclusions. 

• The proposed method is robust to the absence of pure 
pixels; 

• The proposed method provides competitive results in 
terms of aSAM while allowing endmember variability 
to be estimated for each endmember in each pixel; 

• For most simulation scenarios, the abundance MSEs and 
the REs are lower than the MSEs and REs resulting from 
state-of-the-art methods; 

• The proposed method is computationally more expensive 
than existing algorithms. 

We can note that the smoothness penalization on the abun¬ 
dances proves to be particularly appropriate in this experiment. 




0.8 
o 0.6 

CO 

o 

^ 0.4 

CD 

a: 

0.2 


500 1000 1500 2000 

X (nm) 

(c) Endmember 3 

Fig. 4. Endmember estimations obtained on synthetic data in absence of pure 
pixels (cf. Figs. 3 for the abundance estimations). The ssmvBCD/ADMM- 
estimated endmembers (red lines) are given with typical examples of the 
estimated variability (cyan dotted lines). The VGA endmembers are given 
in blue dotted lines for comparison. 



Moreover, an increasing number of endmembers implies a loss 
of estimation performance. This result can be expected since 
VCA/FCLS algorithm is used as an initialization step. 

Finally, the variability captured by the proposed model 
is presented in Figs. 3 and 4 for three endmembers: the 
difference between the variability intensities detected in the 
upper and the lower part of the scene is due to the different 
variability coefficients applied to these areas, thus illustrating 
the consistency of the proposed method. 

V. Experiment with real data 
A. Description of the datasets 

The proposed algorithm has been applied to real hyper- 
spectral datasets obtained by the Airborne Visible Infrared 
Imaging Spectrometer (AVIRIS). The first scene was acquired 
over Moffett Field, CA, in 1997. Water absorption bands were 
removed from the 224 spectral bands, leaving 189 exploitable 
spectral bands. The scene of interest (50 x 50) is partly 
composed of a lake and a coastal area. 

The second scene is a 190 x 250 image extracted from the 
well-known Cuprite dataset^. The number of spectral bands is 
189 after removing the water-absorption and low SNR bands. 
Many works previously conducted on this image provide 
reference abundance estimation maps. 

The parameters used for the proposed approach are identical 
to those used for the experiments with synthetic data (see Table 
II). The only difference is that the algorithm has been stopped 
when the relative difference between two successive values 
of the objective function is less than 10“^. This value has 

^The Moffett and Cuprite images are available at http://www.ehu.es/ 
ccwintco/index.php?title=Hyperspectral_Remote_Sensing_Scenes, and http:// 
aviris.jpl.nasa.gov/ 
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TABLE III 

Simulation results for synthetic data in presence of pure pixels (GMSE(A)x 10“^, GMSE(dM)xlO“^, RE xlO“^, 7 = !)• 


K = 3 , (a, 13) 

= (1.4, 2.5 X 10-5) 

K = e, {a, 13) = (0.37, 5.1 x lO-^) 

aSAM(M) C) GMSE(A) 

GMSE(dM) RE 

time (s) 1 aSAM(M) (°) GMSE(A) GMSE(dM) RE time (s) 


VCA/FCLS 

6.0038 

3.80 

/ 

7.56 

1 

6.3313 

2.24 

/ 

2.92 

1 

SISAL 

5.2665 

3.08 

/ 

3.35 

2 

3.8365 

3.05 

/ 

2.25 

3 

FDNS 

6.0038 

3.79 

/ 

7.56 

4 

6.3313 

2.22 

/ 

2.92 

5 

AEB 

5.6971 

2.07 

/ 

3.50 

52 

5.7017 

1.31 

/ 

2.40 

142 

BCD/ADMM 

5.9910 

3.51 

4.00 

0.20 

92 

6.2965 

1.59 

2.93 

0.05 

230 

ssBCD/ADMM 

5.7765 

3.15 

4.25 

0.23 

422 

6.0304 

1.44 

2.97 

0.07 

848 

s smvB CD/ADMM 

5.4390 

3.01 

4.25 

0.25 

530 

6.3397 

1.42 

2.97 

0.07 

603 


TABLE IV 

Simulation results for synthetic data in absence of pure pixels (GMSE(A)x 10“2, GMSE(dM)xlO“^, RE xlO“^, 7 = 1). 


K = 3, {a, p) = (24.5,4.2 x IQ-^) 

K = &, (a,P) = (0.71,4.8 x IQ-^) 

aSAM(M) (°) GMSE(A) GMSE(dM) RE 

time (s) 1 aSAM(M) (°) GMSE(A) GMSE(dM) RE time (s) 


VCA/ECLS 

5.0639 

2.07 

/ 

2.66 

1 

6.5530 

2.52 

/ 

2.82 

4 

SISAL 

4.4318 

2.16 

/ 

2.56 

2 

6.0431 

1.63 

/ 

2.02 

5 

FDNS 

5.0639 

2.06 

/ 

2.66 

3 

6.5530 

2.53 

/ 

2.82 

7 

AEB 

5.1104 

2.11 

/ 

2.66 

33 

6.0016 

1.78 

/ 

1.85 

208 

BCD/ADMM 

5.2480 

2.13 

3.81 

0.25 

140 

6.2785 

2.14 

3.33 

0.30 

3041 

ssBCD/ADMM 

4.1549 

1.44 

4.36 

0.38 

1263 

6.2763 

1.74 

3.04 

0.076 

1527 

ssmvBCD/ADMM 

5.0584 

1.94 

4.59 

0.47 

1667 

6.3207 

1.67 

3.05 

0.08 

795 


been chosen to obtain a compromise between the estimation 
accuracy and the computational cost implied. The values 
selected by cross-validation for a, P and 7 are given in Table 
V. 

TABLE V 

Experiment results conducted on real data [ssBCD/ADMM for 
Moffett with (a, (3) = (0.05,0), ssvcaBCD/ADMM for Cuprite 
WITH (q;,/3) = (0.014,404), RE xlO -^,7 = 1]. 


Moffett Cuprite 



RE 

time (s) 

RE 

time (s) 

VCA/FCLS 

2.50 

0.4 

3.69 

9.9 

SISAL 

1.12 

30 

2.16 

15 

FDNS 

2.69 

1 

3.69 

11 

AEB 

6.25 

10 

0.40 

615 

BCD/ADMM 

0.18 

144 

0.23 

1.9e4 


B. Results 

The unmixing performance are reported in Table V. For 
the Moffett image, the variability detected by the proposed 
algorithm is displayed in Figs. 5 and 6 . The variability seems 
to be more significant on the coastal area where the mixture 
is not appropriately described by a linear model. The potential 
non-linearities usually observed close to the coastal areas [32]- 

[34] are interpreted as variability in the proposed method, 
which tends to corroborate its consistency. Note that the 
advantage of the proposed method is that it does not require 
to consider a sophisticated non-linear model accounting for 
interactions between the different endmembers as in [32], [34], 

[35] . Conversely, all deviations from the LMM are contained 
in the variability components dMn,/c. We can also note that 


the variability peaks observed in Fig. 6 are a clear indication 
that several corrupted spectral bands have not been removed 
prior to the unmixing process. 

The results obtained for the Cuprite scene are reported in 
Figs. 7, 8 and 9. Comparing our results with those of [28], we 
visually found out that some similar endmembers that were 
identified as different signatures by VCA for K = lA [28] 
are interpreted as multiple instances of single endmembers 
in our setting (K = 10). The identification is given in 
Figs. 7 and 9. Fig. 8 shows that the algorithm captured a 
significant variability level in the pixels where many different 
endmembers are detected, which reveals that the spectral 
mixture may not be strictly linear in these pixels. 

VI. Conclusion and euture work 

This paper introduced a new linear mixing model accounting 
for spatial and spectral endmember variabilities. The proposed 
model extended the classical LMM by including an additive 
spatially varying perturbation matrix that can capture endmem¬ 
ber variabilities. The resulting unmixing problem was solved 
by alternating marginal minimizations of an appropriately 
regularized cost function, each minimization being performed 
by an ADMM algorithm. Simulations conducted on synthetic 
and real data enabled the interest of the proposed solution 
to be appreciated. Indeed, the proposed method compared 
favorably with state-of-the-art approaches while providing a 
relevant variability estimation. The choice of the penalization 
parameters P and 7 was performed by cross validation. 
We think that it would be interesting to develop automatic 
strategies for estimating these parameters. Finally, due to 
the significant number of unknown parameters, the proposed 
method is not intended to be applied to very large images. The 
proposed approach can be applied as a complementary tool 
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when analyzing small hyperspectral images a priori believed 
to be affected by a non-negligible variability level. Decreasing 
the computational complexity of the algorithm introduced in 
this work is clearly an interesting prospect. 


Appendix A 

Constraints and Penalization Terms 

A. Abundance penalization: spatial smoothness 

The abundance smoothness is expressed in matrix form as 


$(A) = i ||AH||2 


(31) 



(a) Abundance I (b) Abundance 2 (c) Abundance 3 


where H denotes the matrix computing the differences be¬ 
tween the abundances of a given pixel and the respective 
abundances of its 4 neighbors 


H = 


I I I H; 




For h = 1,... ,H, we introduce 



(d) Abundance I (e) Abundance 2 (f) Abundance 3 



(g) Variability 1 (h) Variability 2 (i) Variability 3 


Fig. 5. ssmvBCD/ADMM abundance estimations (Moffett scene). 
VCA/FCLS results are shown in Figs. 5a to 5c whereas Figs. 5d to 6c are for 
the ssmvBCD/ADMM. The spatial distribution of the variability with respect 
to each endmember is presented in terms of energy ||dm^^fc ||2 for the 
kth endmember in the nth pixel) for visualization purpose. 


/O -1 0 

0 1 


H 


h — 


0 \ 


0 


])WxVL 


1 -1 

Vo . 0 ij 

/I 0 . 0 \ 

-1 1 


H 


h = 


0 

Vo 


1 0 

0 -1 0 / 




Hence 


=Diag(Hi,...,HH) and = Diag(Hi,..., H^). 


In addition 


= [0Ar,iv,H|] and = [H|,0Ar,w] 


with 



500 1000 1500 2000 2500 

X (nm) 

(a) Endmember 1 




Hj = 


Hi = -Hj. 



^ -1 0 ... 

0 

\ 

w 





1 

0 



0 

-1 


N -W 





^ 0 ... 0 

1 

/ 


The only terms in ^ ||AH||p related to are 




^n,n-\-kN 


k=0 


||2 


X (nmj 
(c) Endmember 3 

Fig. 6. ssmvBCD/ADMM endmembers (Moffett scene). The 

ssmvBCD/ADMM-estimated endmembers (red lines) are plotted with 
the VC A endmembers (blue lines) for comparison, and typical examples of 
the estimated variability are given in cyan dotted lines. 


cAri 
N 3 


EE hn^n+kN^i^n+kN^i j 


i=l k=0 
i^n 


^NxiN-W) 


( 32 ) 
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B. Endmember penalization 

1) Distance between the endmembers and reference signa¬ 
tures: The distance between the endmembers and the available 
reference signatures is 

1 1 ^ 

^-(M) = - ||M - Moll^ = - ||m, - mo,dl 2 • (33) 

^ E=1 

As a consequence, the penalty for the ^th band is 

^ ||m£ - mo /||2 . (34) 

2) Mutual distance between the endmembers: The distance 
between the different endmembers can be expressed as follows 


b) Positivity constraint on M: Using the following nota¬ 
tions 

Y = UYp,oj.+Yi, Yi = [y|...|y] 

M = UT + Y 2 , Y 2 = [y|...|y] gR^^^ 

one has 

'^ir — ^ ^ n^jtjrp yi — ^ ^ uijtjrp -j- -|- y£. 

3 

The positivity constraint for mir can then be expressed as 

^ Vi 

tkr ^ • 


K / K 


*Jll 2 


= Hz ( 

i=i Vi=i / 

= if;i|MGfe||2 = i||MG| 


(35) 


k=l 


Introducing the two sets of integers 

^k = > 0 } 

Uf^ = {i\u£k < 0 } 

the previous equation implies that tkr ^ 


with 


G = 


Gi 


G 


and for /c = 1 ,..., if 

G/c = —Ik 


K 


■ekll 


:)KxK^ 


hr = max - 

eeu+ \ 


= min I - 


/cr’ ^kri 

yi + '^j^k 
U£k 

yi + 


£eu~ \ U£k 


(38) 

(39) 


where ek denotes the kth vector in the canonical basis of 
Hence 

K 

|2 


V’(m0 = t'^Wm.tQkWl. 


iK 


(36) 


4) Positivity constraint on M and dM: This case differs 
from the previous one as the positivity constraint must be 
verified simultaneously by M and = M + dM^. We 
will consequently derive a condition similar to (38). Let 
be the projection of in the PC A subspace 


k=l 


3) Volume and endmember positivity constraint: The vol¬ 
ume penalization is expressed using T, hence the need to find a 
condition equivalent to the positivity of M (see [36]). We will 
first analyze the general expression of the volume penalization 
with respect to t/., and then give a condition on T ensuring 
the positivity of M (respectively M + dM^ when endmember 
variability is considered). 

a) Volume: The determinant of a matrix X G 
can be developed along its Ah row yielding 

det(X) = det(Xij) = 


M„, = UT„ 


Since 


with 


T„ = V (M„ - Y 2 ) = (T + dT„) + VY 2 


T = V (M - Y 2 )_ 
dT„ =V(dM„-Y 2 ) 

the positivity constraint can be written 


with 


m^^>0 4^ hr>-- 


U£k 


f,= [(-ir+^det(Xy)].^^ G 
Consequently, for /c = 1,..., if — 1 




with 


det 




= t/ef/c- 


Using previous developments 
i’hk) = 


1 


2{K-iy? 


(tfcffc)^ 


(37) 


^ hr ^ [hr hit] 


^kr ^kr 4 “ ^^kr 4 ~ ^kr 

^ ye + T^j^k£hi3: 

m + 'Zj^k'^^ihr 


.n- 

^kr 


= max I — - 

ieu+ 


^n+ 

^kr 


mm 

£eU7: 


U£k 


(40) 

(41) 

(42) 
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We introduce the functions gk defined by 

^Ix/sT _ y |^2(N+l)x/sr 


9k ■ 


-^+^k 


ljv(x + Zfc) + 


-ljv(x + Zfe) + 


^ dtl 

— tjv,jk/ 
t+j, - dti,fe ^ 

VJV.fc “ ^^N,k) _ 


where 


t+ — [/■+ . 


’ ^kK\ 

’ ^kx]' 


^(A) (^a„, i ||y„ - (M + d]V[„)a„ 

+ I {cAn ||a„||2 + 2cja„) +1^+ 


+ 


Mn 


(A) 


Qa„ + - s + 


B. Resolution with respect to M 

1) Distance between the endmembers and reference spec¬ 
tral signatures: Using (34), the scaled augmented Lagrangian 
(26) is 


U'i 


(m<!, 


(M) .(M)\ 








/3 


em^ - 


ye - A - de 

2 
F 


,(M) 


+ -||m^-m^,o|l2 


^s+ 


(wr>) ^ 


Thus 


m, = 


(43) 


- Se^ 

{' 


■ I3me,0+ 




(49) 


AA^ + ^r^ (e^e + /?) I 


1 -1 


K 


and for k = 1,.... K 


Finally the positivity constraint on the sum of the endmembers 
and their variability can be written 

nil + dm^ ^ 0 ^ V^, Vn (44) 

^ gk{tk) h 02(iv+i),K VA: = 1 ,..., iT - 1 . (45) 

C Variability penalization 

The variability energy penalty is 

T (dM) = i ||dM||2 ^ V (dM„) = i ||dM„||2 . (46) 

Appendix B 

Solutions to the optimization sub-problems 
A. Resolution with respect to A 

Using (32), the scaled augmented Lagrangian (23) becomes 


w 


(M)* 

e,k 


= max I [em^ 


.(M)i 


c, Oat+I^ . 


(50) 


In the absence of any endmember penalization, the solution is 
obtained by making /3 = 0 in the previous equation. 

2) Mutual distance between the endmembers: Using (36), 
the scaled augmented Lagrangian (26) is 


(M) 




W 


(M) a(M))^1 

(M) 

~ ’nsr(lVI) , 

em<! - 4 


ye - m^A - 5e 
2 
F 


Ar 


K 


^Y.\\ikeGu\\l+Is^ (W^^). 


k=l 


Thus 

- 5l< 

nil 


IW 


A" + - F, - 


K 


AA^ +/3J2 GkGl + 4^^ (e^e) I 


£- ^^l J 
-1 
K 


k=l 


(51) 


with given by (50). 

3) Volume penalization: Since the penalty is expressed 
with respect to the variable T, the optimization sub-problems 
related to the endmembers have to be re-written accordingly. 
-1 Using the notations 

At = [ dTiai | ... | dT^raAr ] (52) 

Yi = [y|---|y] (53) 


= 


i-K 


Thus, for n = 1,..., 

(M + dM„)^(M + dM„) + + «cA„I 

(M + dM„) V - «c„ + (s - - Ai^)) 

(47) 

and 

= max Oif j (48) 

where is the vector composed of the K first elements of 

and the max must be understood as a term-wise operator. 
In the absence of any penalization, the solution is obtained by 
making a = 0 in the previous equations. 


Y 2 = [y|-"|y] e 


^LxK 


(54) 


we obtain 

||Y-MA-A||2= ||U(Yp,oj. ... 

= ||U(Yp,oj.-TA-AT"' 


TA-At)+Yi - 2 Y 2 AI 

^r)iiF 


2(U(Yp™j.-TA-A t) 


Yi - 2 Y 2 A 


Yi - 2 Y 2 A 


F ■ 
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The only terms depending on T are 

||Yp,oj, - TA\\l + 2 /At- U^(Yi - 2Y2A) 


TA 


with 

N .K-1 . 

{S|TA) = Tr (S^TA) = ^ ( E ) ' 

n=l^j = l / 

For /c = 1,..., A — 1, the resulting sub-problems are 


tl = arg min < 


y^' - t/cA ^ + T.n=li^kntk^n) 

2(K-1)!2 i^k^k) 

S.t. gki^k) E 02(Ar+l),K 


> . 


fT) 

Introduce the splitting variables ^ such that 


Vfc = l, 


,K-1. 


According to (37), the scaled augmented Lagrangian is 

£,m(t,,wf\Af)) = 




+ 


N 


/3 


n=l 

2 


2(A-1) 


T2(tA) +x5+(wr) 


+ 


(T) 


fffc(tfe)-W["' + A 


(T) , . (T) 


Finally, 


t* — 




n=l 


AT 


U/c + + X](^n,/c + '^n,/c) 


n=l 


+/.r[i -1 -isi (wr> - Ar>) 


AA^ + 




-| -1 


(A - 1)!2 

where Z = VY 2 and for p = 1,..., A 

( 5fe(tfc) + , 02(JV+1) ) • 


C Resolution with respect to dM 

Using (46), the scaled augmented Lagrangian (28) is 


Hence 
dM* = 


(y„ - Ma„)a^ + (w<,<I“> - M - AW^I) 


ana^ + 


and for p = 1 ,..., A 

= max (^[dMn + M + 0^) • 


(59) 



(a) Sphene 


(b) Alunite 


(c) Dumortierite 


(55) 



(d) Montmorillonite (e) Andradite 


(f) Pyrope 


(56) 



(g) Buddingtonite 


(h) Muscovite 


(i) Nontronite 



(57) 


(j) Kaolinite 

Fig. 7. ssvcaBCD/ADMM abundance results (Cuprite scene). The given 
identification is based on a visual comparison with the results obtained in 
[28]. 



(58) 


(a) Sphene (b) Andradite (c) Muscovite 

Fig. 8. Spatial distribution of the variability with respect to the endmembers 
presenting the most significant level of variability. The variability is presented 
in terms of energy (-^ ||dm^ ||2 for the Aith endmember in the nth pixel) 
for visualization purpose. 


>C^(dM) AdMn^ — 
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